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Abstract 

In this paper we study the motion of three linked ellipses moving through a 
viscous fluid in two dimensions. The angles between the ellipses change with time in 
a specified manner (the gait) and the resulting time varying configuration is similar 
to the appearance of a swimming leech. We simulate the motion using the particle 
method Smoothed Particle Hydrodynamics (SPH) which we test by convergence 
studies and by comparison with the inviscid results of Kanso et al. (2005) and the 
viscous results of Eldredge (2006, 2007, 2008). We determine how the average speed 
and power output depends on the amplitude and oscillation frequency of the gait. 
We find that the results fit simple scaling rules which can related to the analytical 
results of G.I. Taylor for the swimming of long narrow animals (1952). We apply 
our results to estimate the speed of a swimming leech with reasonable accuracy, 
and we determine the minimum power required to propel the bodies at a specified 
average speed. 



1 Introduction 



The subject of this paper is the motion of linked rigid bodies moving in a weakly com- 
pressible, viscous fluid. It is closely connected with mathematical and computational 
studies of the swimming of fish, and with the motion of underwater vehicles and robotic 
fish propelled by changes of shape. Our approach is primarily computational using the 
SPH algorithms of Kajtar and Monaghan (2008) to establish scaling relations for the 
motion. The work is closely related to recent work on the motion of linked bodies in an 
infinite, two dimensional fluid which may be inviscid (Kanso et al., 2005, Melli et al., 
2006) or viscous (Eldredge, 2006, 2007, 2008). When the fluid is inviscid it is possible to 
bring powerful mathematical formalisms to bear on the problem in a manner similar to 
the motion of a single body in an inviscid fluid (see for example Lamb, 1932). However, 
for problems involving free surfaces, or complicated rigid boundaries, or a stratified fluid, 
these methods become very complicated. Our approach is capable of handling arbitrary 
body shapes and boundaries though in the present paper we concentrate on motion of 
linked ellipses moving in a periodic domain. The SPH method also has advantages over 
the vortex particle method of Eldredge (2006) for problems where the bodies penetrate 
a free surface, but in the present case no such difficulties exist and the vortex particle 
method provides a convenient comparison for the SPH calculations. 

The bodies we consider are solid bodies linked by virtual rods which join at pivot 
points. The rods are described as virtual because they do not have any mechanical 
function except to define the direction of fixed lines in the bodies. In particular, fluid 
can flow between the ellipses. The angles between the rods (and therefore the bodies) are 
specified as an oscillating function of time. A specification of the time variation of these 
angles is called the gait. 

Our aim is to determine the scaling relations which relate the speed and power 
developed by the linked bodies to the frequency and amplitude of a standard gait which 
propels the linked bodies along a path which is, on average, a straight line. A related 
problem was considered by Taylor (1952) who studied the motion of a long slender body 
and applied his results to the motion of a leech and a snake moving in water. Our 
three ellipse system has a motion which is similar to that of the leech and snake because 
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their oscillations are are roughly sinusoidal, and are therefore not too different from the 
oscillations of our connected ellipses. Taylor's analysis provides a remarkably accurate 
guide for the functional form of the dependence of the velocity and power of our three 
ellipse system on the frequency and amplitude of the gait. 

The plan of the paper is to first discuss the SPH algorithm. We then show by 
convergence studies that a periodic domain can be used to represent the infinite domain 
with errors that are typically 5%. We also establish convergence with resolution. We first 
compare our results with the viscous results of Eldredge (2008) for massless bodies by 
a series of simulations with decreasing body mass. The agreement is very satisfactory. 
We then compare our results to the inviscid results of Kanso et al. 2005 by changing the 
viscosity so that the Reynolds number varies from 50 to 5000. This comparison shows 
that the SPH results converge to the inviscid results for the highest Reynolds number. 
We then discuss scaling relations for the velocity and power output, and relate them 
to Taylor's (1952) analytical relations for the velocity and power output of long narrow 
animals swimming. We apply these results to the swimming of a leech. Finally we 
determine the minimum power, and corresponding gait, required to propel the bodies at 
a specified average speed. Throughout this paper we use SI units. 

2 Equations of motion and constraints 
2.1 Equations of motion 

We consider motion in two dimensions and use cartesian coordinates. A typical config- 
uration of the bodies is shown in Figure [T} The motion of the fluid, which is assumed 
incompressible, is specified by the Navier Stokes equations. In cartesian tensor form these 
equations are 




(2.1) 



where Nb denotes the number of bodies, is the stress tensor 
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P is the pressure and /x the shear viscosity coefficient. The function 5(sk) is a one di- 
mensional delta function, and Sk is the perpendicular distance from the surface of 
body k to the position where the fluid acceleration is required. The unit vector rij(k) 
is directed from body k into the fluid. The introduction of forces into the acceleration 
equation ( 2.1| as an alternative to specifying boundary conditions on the velocity is due 
to Sirovich (1967, 1968). In his formulation, as in ours, S(s) is a delta function defined so 
that for any quantity B(r) 

r BS(s)dv = [ BdA, (2.3) 



where the first integral is over the volume and the second integral is over the surface. 
In this way a volume integral involving a delta function becomes equivalent to a surface 



integral, and the body force per unit volume in (2.1) becomes a force per unit area. 
This force provides both the pressure which prevents penetration of the rigid body, and 
the viscous traction term. It mimics the fundamental molecular basis of the boundary 
conditions namely that the atoms of the fluid do not penetrate the atoms of the solid 
because of the atomic forces between the liquid and the solid atoms. 

A closely related method of using boundary forces is due to Peskin (1977) who sim- 
ulated elastic membranes such as the heart interacting with a fluid. Peskin's equations 
(2.3) to (2.6) are essentially those of Sirovich, though the Peskin deals with an elastic ma- 
terial and Sirovich assumes the body is rigid. Further details about Peskin's formulation 
can be found in Peskin (2002). 

We denote an element of area on the surface of body k by dA(k). The motion of 
the centre of mass R(/c) of solid body k (with mass M(k)) is given by 

M{k) jf± =- a ij n 3 {k)dA{k) + F\k), (2.4) 

where F(k) is the force due to the constraints. The rotation of rigid body k, with moment 
of inertia I(k)), is given by 

m ^ = / (dfc X h)dA{k) + T(fc) ' (2 - 5) 

where d& is a vector from the centre of mass of body k to the element of area dA^, b is 
the force on the element of area, and r(k) is the constraint torque on body k. 
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Figure 1: The configuration of the bodies (assumed to be ellipses). The link position is denoted by a 
filled circle. The straight line through a body passes through its centre of mass and is assumed to be 
rigidly attached to the body with one end attached to the link. The angles 8 are defined relative to a 
fixed direction in space (shown by parallel dotted lines) which is taken to be the x axis of a cartesian 
coordinate system in our calculations. The angles ip determine the gait and are specified functions of 
time. 

In the following, to simplify the notation, the subscript k will always denote the 
label of a body. Thus, for example, H(k) will be replaced by 

2.2 Constraint equations forces and torques 

The angle 8k which fixes the rotation of body k is defined as the positive rotation of a 
line fixed in the body from the x axis of a cartesian coordinate system fixed in space. For 
simplicity we assume the line fixed in the body is an axis of symmetry. The constraint 
conditions on the angles are 

V?m = Orn+l ~ m , (2.6) 

where m is the link number and (p m is a specified function. The form of the (p m determines 
the gait of the bodies. For the examples we consider here there are three bodies and two 
links as shown in Figure [T} In the simplest case ip m is a function of t but, in general, it 
depends on other variables. For example, in a biological problem, it could depend on the 
centre of mass coordinates in such a way that the fish slows down when it enters a region 
where food is abundant. 

In addition to the constraints on the angles there are constraints associated with the 
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links. We assume the link, or pivot, is at a distance 4 from the centre of mass of body k. 
The condition on the X components of the centres of mass of bodies k and k + 1 is that 
the X coordinate of the link between them is given by 

X k - 4 cos (0 fc ) = X k+l + 4 +1 cos (0fc+i), (2.7) 

or 

X fe - 4 cos (6> fe ) - X fe+ i - 4+i cos (0 fc+ i) = 0. (2.8) 
Similarly the Y constraint is 

Y k - 4 sin (6 k ) - Y k+1 - 4+i sin (6 k+l ) = 0. (2.9) 

These constraints enable the coordinates of the centres of mass of the bodies, and their 
angles 9 to be written in terms of those of any selected body. Similarly, by differentiating 
the constraint conditions with respect to time, the velocities X and Y and angular velocity 
f2 of the bodies can be written as functions of the same selected body. The number of 
degrees of freedom (coordinates and velocities) of N linked bodies in two dimensions is 
therefore 6 compared with the QN degrees of freedom of N independent bodies in two 
dimensions. If the ip m are functions of t alone it is possible to reduce the equations of 
motion to those involving the coordinates and velocities of one of the bodies. This can also 
be done when the <p m are functions of both coordinates and time but it is inconvenient to 
eliminate variables and, in our view, simpler to take account of the constraints by using 
Lagrange multipliers. For that reason we use Lagrange multipliers even though, in the 
applications to be described in this paper, the ip m are functions of t only. 

For the case of three bodies we have two links and therefore 6 constraints. We 
denote the Lagrange multipliers for the X, Y and 9 constraints of link m by A^, Ay^ 
and Ag"^ respectively. Using standard methods for holonomic constraints (e.g. Landau 
and Lifshitz, 1976) we find the following expressions for the constraint forces and 
torques r k for the various bodies. For body 1 

F 1 = (A«,A^ 1) ), (2.10) 

for body 2 

F 2 = (-A?,-A?) + (Ag ) ,Ag ) ), (2.11) 
7 



and for body 3 

F 3 = (-A?,-Ag ) ). (2.12) 

These constraint forces do not affect the total linear momentum of the bodies because 
they sum to zero. 

The constraint torque on body 1 is 

Tl = -\V> + A£^i sin (0i) - Ag^i cos {9 1 ), (2.13) 

on body 2 it is 

r 2 = A< 1} - A< 2) + (a^ } + A<?) £ 2 sin (0 2 ) - (a^ 1} + A^) £ 2 cos (0 2 ) (2.14) 
and on body 3 it is 

r 3 = A^ + A?£ 3 sin (0 3 ) - A? } £ 3 cos (0 3 ). (2.15) 

The constraint forces and torques are provided by the engines which drive the an- 
gular variation between the bodies. In the case of fish these engines are the muscles of 
their bodies and the work done is provided by the internal chemical energy generated by 
the fish. The way these constraint forces and torques affect the angular momentum will 
be discussed in §5. 



3 SPH equations for the fluid 

The form of the SPH equations that we use is discussed in more detail by Monaghan 
(1992, 2005). For the liquid SPH particles the acceleration equation is 

dv f P P \ Nb 

-IT = - E mb NT + "4 + n «6 VaWab + EE ^ " mjUajVaWaj] . (3.1) 
at b \Pa Pb J k=ijes k 

In this equation the mass, position, velocity, density, and pressure of particle a are m a , 
r a , v a , p a , and P a respectively. W ao denotes the smoothing kernel W(r a — r b , h ab ) and V a 
denotes the gradient taken with respect to the coordinates of particle a. In this paper W is 
the fourth degree Wendland kernel for two dimensions (Wendland, 1995), and has support 
2h ab . In the present calculations the h ab used in W ab is an average h ab = (h a + h b )/2. The 
choice of h is discussed in detail by Monaghan (1992, 2005). In this paper we choose h 
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to be 1.5 times the initial particle spacing so that the interaction between any two fluid 
particles is zero beyond 3 initial particle spacings. 



The first summation in (3.1) is over all fluid particles and is the SPH equivalent of 



the first term on the right hand side of (2.1). The last term in (3.1) is the contribution 
to the force per unit mass on fluid particle a due to boundary particles and is equivalent 



to the last term in (2.1). A body label is denoted by k, and j G Sk is one of the set of 
boundary particle labels on body k. The term f OJ - is the non-viscous boundary particle 
force per unit mass on fluid particle a due to boundary particle j. In the present paper 
we use the boundary forces analysed by Monaghan and Kajtar (2009). The force f a j acts 
on the line joining particle a and j. The boundary particles delineate the boundaries, and 
produce forces on the fluid in a similar manner to the delta function forces of Sirovich 



discussed after (2.2). 

The viscosity is determined by 11^ for which we choose the form (Monaghan 1997, 

2005) 

Pa&| r a&| 

In this expression a is a constant, and the notation v a b = v a — is used. p ab denotes the 
average density \{p a + Pb)- We take the signal velocity to be 

V sig = \{c a + C b ) - 2 V " b Tab , (3.3) 

2 r ab 

where c a is the speed of sound at particle a (Monaghan, 1997, though here we take v sig 
to be half used in that paper and a is therefore a factor 2 larger). v s i g is dominated by 
the terms involving the speed of sound. The kinematic viscosity can be estimated by 
taking the continuum limit which is equivalent to letting the number of particles go to 
infinity while keeping the resolution length h constant. By a calculation similar to that 
in Monaghan (2005) it is found that the kinematic viscosity for the Wendland kernel is 
given by 

v = -ahv si g. (3.4) 

o 

SPH calculations for shear flow agree very closely with theoretical results using this kine- 
matic viscosity (Monaghan, 2006). 
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The pressure is given by 



2 / / \ 7 

Po C a I Pa 



K ( [fj -1). (3.r,) 

where p is the reference density of the fluid. To ensure the flow has a sufficiently low 
Mach number to approximate a constant density fluid accurately, we determine the speed 
of sound by c a ~ 10 V where V is the maximum speed of the fluid relative to the bodies. In 
this dominated by the first two terms. The precise value of c a will be specified 

for each simulation. 

The form of the SPH continuity equation we use here is 

^ = $> b v a6 ■ VWU (3.6) 

b 

and the position of any fluid particle a is found by integrating 

dF = v " (3 ' 7) 

In the present simulations the liquid SPH particles were initially placed on a grid 
of squares and thereafter allowed to move in response to the forces. The time stepping of 
the SPH equations uses an algorithm which is symplectic in the absence of dissipation. 
The details of this scheme are given by Kajtar and Monaghan (2008). 

4 SPH equations for the rigid bodies 

The non-viscous force on boundary particle j due to all fluid particles is 

f^=m,£f, a , (4.1) 

a 

where f} a is the force per unit mass on boundary particle j due to fluid particle a. The 
viscous force is 

= -m-j m a U aj VjW aj = rrij ^ m a U aj W a W aj , (4.2) 

a a 

where we have used the fact that VjW a j = —V a W a j. The total force on particle j is 

f. = f M +f M. (4.3) 
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The equation for the centre of mass motion of body k is then 

iV 



jes k 



and the torque equation is 

dt 



4 ^ = ^ (r ._ Rfc)xf . + Tfe (45) 



The motion of a boundary particle can be determined from the motion of centre of 
mass and the rotation about the centre of mass. Thus for particle j on body k, 

dr- 

^=V t + fl fc zx ( Tj -R fc ), (4.6) 

where, in this two dimensional problem, the rotation is around the z axis which is per- 
pendicular to the plane of the motion. 



5 Conservation of linear and angular momentum 

The total rate of change of the linear momentum of the rigid bodies with respect to time 
is 

E M ^ = E E ^ = E E E m > ft- - m a v 3 w a] \ (5.i) 

k k jeSk k jGSk a 

where, as noted earlier, the sum over the constraint forces is zero. The rate of change of 
linear momentum of the fluid SPH particles is given by 

E ma ~dT = E E E ma [ fa J ~ rrijUajV a W aj \ . (5.2) 

a a k j&Sk 

noting that the sum over the pressure and viscous forces between fluid particles vanishes 
because of symmetry. 
Recalling that 

^2 m &a = ~y] mj aj and V a W aj = -VjW aj , (5.3) 

a a 

we deduce that 

V M,' . 

dt ^ dt 



k 
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which shows that the linear momentum 

(5.5) 

k a 

is conserved. 

The angular momentum of the bodies is composed of the centre of mass angular 
momentum about some fixed origin, and the sum over each body of the spin angular 
momentum about the centre of mass of body. The time rate of change of the total centre 
of mass angular momentum is 

k k j£S k k 

The rate of change of the spin angular momentum is 

E^ = EE^- R *) x f i + E r *- ( 5 - 7 ) 

k k j£S k k 



The rate of change of the angular momentum of the fluid particles is 

dv, 



£ m a r a x ^ = £ £ £ maVa x [ f< w ~ mjIlajW a W aj ] , (5i 

a k j£S k 



where, because of symmetry, the sum over pressure, and viscous terms between fluid 
particles have vanished. The rate of change of the total angular momentum (the sum of 



(5.6), (5.7) and (5.8)) becomes 

a k jeS fc 

+A^ (-Fi + i x sin {fix) + Y 2 + £ 2 sin (0 2 )) 
+A^ (Xi - £ x cos {fix) -X 2 -£ 2 cos (0 2 )) 
+A^ (-y 2 + £ 2 sin (6 2 ) +Y 3 + £ 3 sin (0 3 )) 

(X 2 - £ 2 cos (0 2 ) +X 3 -£ 3 cos (0 3 )) . (5.9) 

The first term vanishes because the boundary forces are radial and the last four terms 



vanish because of the constraint conditions (2.8) and (2.9). 



Finally we note that the previous arguments about conservation assume the time 
derivatives are exact. The actual conservation in the numerical simulations depends 
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on the form of the time stepping algorithm. Linear momentum is always conserved to 
round off error, but the angular momentum conservation is less accurate because the 
Lagrange multipliers are calculated at the mid-point. In our simulations we use periodic 
boundaries and these do not conserve angular momentum exactly. A detailed discussion 
of the conservation of angular momentum is given by Kajtar and Monaghan (2008). 

5.1 Remarks concerning external boundaries 

The SPH equations can be applied to the linked bodies moving in a channel, as is the 
case for many laboratory experiments on fish, or in a pond with an irregular boundary, 
by replacing the boundaries of the pond by boundary force particles as we have done 
for the rigid bodies. The SPH algorithm does not need to be changed if the linked 
bodies move through and out of a free surface, which would be required to mimic the 
motion of dolphins. This facility was used earlier for bodies hitting the water (Monaghan 
and Kos, 2000, Monaghan et al., 2003). In the present paper, where we compare our 
results with those of Kanso et al. (2005) and Eldredge (2008), we need to deal with 
an infinite medium. This cannot be done directly because it would require infinitely 
many particles. One alternative, and the simplest, is to replace fluids of infinite extent 
by periodic boundary conditions. These boundaries alter the solutions of the differential 
equations but the effects are small if the periodic cells are sufficiently large. We determine 
their effect by carrying out test calculations for successively larger domains. 



where throughout this section, we set (3 — 1 and u — 1. The ellipses have semi- major 
axis a = 0.25, semi-minor axis b = 0.1a, and distance between the tip of the ellipse 
and the pivot c = 0.2a. These dimensions, and the gait, are identical to those of Kanso 



6 The motion of the linked bodies 



We consider ellipses moving with the gait 



<Pi = 2 (O)-0i(O)+/3(cos(wt)-l), 
<P2 = 9 3 (p)-e 2 (0)+Pam(wt), 



(6.1) 
(6.2) 
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Figure 2: The configuration of the bodies at time intervals separated by 27r/3 with time increasing from 
top to bottom. 
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et al. (2005) and Eldredge (2008) but we use a different notation for the angles. The 
configuration of the ellipses is shown in Figure [2] at intervals of 1/3 of a period. 

We define the Reynolds number 3? by using the characteristic velocity V = 2au and 
the characteristic length scale L = 2a, so that 

K=^. (6.3) 

The speed of sound c s = 20au, and the boundaries of the ellipses were defined by boundary 
particles with spacing dp/4 (Monaghan and Kajtar, 2009). The motion takes place in a 
domain with periodic rectangular cells. 

The motion of the linked bodies is characterised by the path followed by the centre of 



mass of the middle body. This path will be referred to as the 'stride path'. The gait (6.1 ) 



and (6.2) is oscillatory with period P = 2tt/uj so that the stride path is oscillatory and 
has the shape of a zig zag. We refer to the straight line distance between two consecutive 
lower points of this zig zag, travelled in time P, as a 'stride length'. The results to follow 
show that the stride length is, in general, not constant, in agreement with the results of 
Eldredge. 

Kajtar and Monaghan (2008) showed that the SPH algorithm gave results in good 
agreement with experiments for a driven oscillating cylinder, and for cylinders freely 
oscillating in a channel flow. In this paper we describe three levels of further tests. The 
first of these is concerned with the convergence as the resolution is refined with a fixed 
periodic cell size, and convergence as the size of the cell is increased with fixed resolution 
(§6.1 and §6.2). The latter is to ensure that our comparison with the results of Kanso et 
al. (2005) and Eldredge (2008) is legitimate. The second level of tests is concerned with 
comparisons with the results of Eldredge by studying the stride length when the mass 
of the bodies is reduced and Kanso et al. by studying the stride length change as the 
viscosity coefficient is increased (§6.3 and §6.4). The third level of tests shows that the 
numerical simulations agree with general scaling relations (§7). 

6.1 Test of the convergence with resolution 

Throughout this section, we use a rectangular domain with periodic boundaries aligned 
with the x and y axes of a cartesian coordinate system. The ratio of the lengths of the 
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sides of the domain, the aspect ratio, is 4:3. The fluid spans from x m - m = to x max along 
the horizontal axis, and from y min = to y max in the vertical axis. The initial coordinates 
of the centre of mass of the middle body were always (X 2 , Y 2 ) = (0.4x max , 0.6y max ). For the 
SPH simulations the periodic boundaries were implemented by copying rows and columns 
of fluid particles Ah in width to the opposite boundary, top to bottom, left to right, and 
vice versa. This process guarantees that the fluid particles of interest in the rectangular 
domain get the correct rates of change in each time step. 

We ran the calculations for initial particle spacing dp = 1/30, 1/40, 1/50 and 1/60. 
The domain was of size x max = 4 and ?/ max = 3. For these tests, the bodies were neutrally 
buoyant and 9ft = 200. The simulations for each resolution were run for the same time. 

The stride paths for the four values of dp are plotted in Figure [3] Note that the 
strides for the lowest resolution (dp = 1/30) are significantly longer than for the other 
three finer resolutions. The paths for dp = 1/40, 1/50 and 1/60 lie almost on top of 
one another although, because of the slight differences in average velocity, the differences 
increase with time and we note that the convergence is not monotonic i.e. the results for 
dp = 1/40 are closer to those for 1/60 than are those for 1/50. However, for the three 
smallest resolutions the relative difference between a stride length of one resolution and 
another is at most 5% (for the third stride). Figure [3] also shows that the direction of the 
path is not sensitive to the resolution. The results of this numerical test indicate that a 
fluid particle resolution of dp = 1/40 is sufficiently accurate to determine the stride path 
in length and direction to within 5%. 

6.2 Test of convergence with domain size 

In order to determine a fluid domain size that adequately represents an infinite domain, 
the calculation was run for a number of periodic cell sizes with fixed dp. We ran the 
calculation with four different domain sizes with the same aspect ratio, x max x y max = 4x3, 
5 x 3.75, 6 x 4.5 and 7 x 5.25. Again, the bodies were neutrally buoyant and 9ft = 200. The 
calculations were run for dp = 1/40, but we found that stride paths varied substantially 
from one case to the next. However, with a resolution of dp = 1/60, the stride paths show 
a smoother trend with increased domain size. 
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The stride paths for the four domain sizes are plotted in Figure |4j Note that the 
strides for the smallest domain are significantly longer than for the other three larger 
domains. These results indicate that a domain of size 5 x 3.75 is close to being sufficiently 
large for modelling an infinite domain. The distances travelled in three strides for the 
different domain sizes, and for the two resolutions are given in Table [TJ These values 
demonstrate the large variation for different domain sizes with dp = 1/40. For dp = 
1/60, neglecting the smallest domain, the maximum relative difference is 3% (between 
the 5 x 3.75 and 7 x 5.25 domains). For dp = 1/40 on the other hand, the maximum 
relative difference is 9% (between the 5 x 3.75 and 6 x 4.5 domains). 




0.5 1 1.5 

X 

Figure 4: The stride paths for different domain sizes with fixed dp = 1/60. The crosses denote the stride 
path for domain size x max x y max = 4x3. Open circles, filled circles and the solid line are for 5 x 3.75, 
6 x 4.5 and 7 x 5.25 respectively. Note that for the purposes of comparing the paths on this plot, the 
stride paths have been shifted so that they all begin at (0,0). 

Based on the numerical tests for the resolution and the domain size, we chose dp = 
1/60 and a domain of size 6 x 4.5 for our subsequent production runs. 
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domain 


dp = 1/40 


dp = 1/60 


4x3 


0.9844 


1.0031 


5 x 3.75 


0.9235 


0.9214 


6 x 4.5 


0.8368 


0.9149 


7 x 5.25 


0.8859 


0.8932 



Table 1: Distance travelled in three strides with different domain sizes, and for two different fluid reso- 
lutions. 

6.3 Comparison with Eldredge 

Eldredge (2008) considers massless bodies, which are inconvenient to use with our algo- 
rithm (the expressions for the body velocities with M\~ = are singular). We can, 
however, observe the trend in the motion of the linked bodies as their mass M — > 0. For 
neutrally buoyant elliptical bodies, the masses are Mq = ponab, and moment of inertia 
/ = Mo(a 2 + 6 2 )/4. We ran the simulation for a number of body masses in the range 
0.5M < M < 5.0M in order to determine a relationship between the mass and the 
stride length. For these simulations dp = 1/60. The Reynolds number 3? = 200 is the 
same as in the calculation of Eldredge. 

Eldredge reports that the massless linked bodies have a stride length of 1.45a. The 
stride lengths from the SPH simulations, as well as the result of Eldredge are plotted in 
Figure [5} The line of best-fit shows that there is a linear trend toward Eldredge's M = 
result. 

The results of Eldredge show that the stride lengths vary from stride to stride. The 
second stride is longer than the first, and the third is longer than the second. Our results 
show a similar behaviour. We find that the second stride length is larger than the first 
typically by ~ 18 — 28%, and the third stride length is larger than the second by ~ 5 — 16%. 
The equivalent results of Eldredge are 20%, and 10% which is similar to our results. For 
the inviscid case (discussed below), Kanso et al. show that the stride length is constant. 

Eldredge estimated the stride length by taking the average of the second and third 
strides, and we followed the same procedure in generating the results of Figure [5j 
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Q g I i i i i I i i i i I i i i i I i i i i I i i i i I i i i i 

1 2 3 4 5 6 

M/M Q 

Figure 5: Stride lengths S, scaled with the length parameter a, for the motion in fluid with constant 
viscosity and fixed gait but a range of body masses. The crosses denote the SPH results and the filled 
circle denotes the result of Eldrcdgc. The line of best-fit is also shown. 
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The vorticity generated by the motion of the linked bodies after time t ~ 25 is 
shown in Figure [6] This plot can be compared to the last frame of Figure 10 of Eldredge 
(2008). Note however, that for this figure M = 0.5M , whereas Eldredge has massless 
bodies. The contours of Eldredge are much smoother than those shown in Figure [6j but 
the main features are recognisable. There is one large, and one smaller eddy near the 
rear of the linked bodies, which are in the same positions as with Eldredge, and there is 
an intense eddy generated by the front body. The stride path is in good agreement with 
Eldredge. 




Figure 6: Vorticity field generated by the swimming linked bodies with M = 0.5Mo. This plot is at time 
t ~ 25. The vorticity contours have values in the range -5 to 5 with 40 levels. The stride path is also 
shown. 



Finally, we note that the vortex particle spacing in Eldredge's simulation is typically 
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a/50 compared to our a/15. Eldredge has 280 panels on each body which is close to the 
244 boundary particles per body in the SPH calculation. 

6.4 Comparison with Kanso 

Kanso et al. (2005) consider the motion of three linked ellipses in an inviscid fluid. 
Although the SPH algorithm is only stable with non-zero viscosity we can study the 
stride length variation with change in the viscosity and estimate the stride length for zero 
kinematic viscosity coefficient as v — > 0. We ran the simulation with neutrally buoyant 
bodies for viscosity in the range 5 x 10~ 5 < v < 5 x 10 -3 , which corresponds to a range 
in Reynolds number of 50 < 3ft < 5000. The calculations were run with neutrally buoyant 
bodies. 

Kanso et al. find that the stride length for the neutrally buoyant bodies is 3.27a. 
This result, as well as the stride lengths from the SPH simulations, are plotted in Figure 
[7| The SPH results show a trend toward the v = case of Kanso. In some respects 
the agreement is remarkable because there are significant differences between the inviscid 
and non-inviscid cases. For example the flow produced by an oscillating cylinder changes 
dramatically as the Reynolds number changes from small 3? ~ 10 to large 3? ~ 1000 
though the time averaged drag terms are nearly constant for 100 < 3? < 10000. In 
the inviscid case the fluid motion produced by a system of oscillating linked bodies will 
cease the instant they stop oscillating, while in the viscous case, the motion will continue 
though it will be damped. And as discussed in the previous section, the strides increase 
in length both for our calculations and those of Eldredge whereas those of Kanso et al. 
are constant. These results suggest that when 3? > 1000 the average motion of the linked 
bodies is determined primarily by added mass effects as discussed by Saffman (1967) for 
swimming by shape change in an inviscid fluid. 

The stride lengths plotted in Figure [7] were calculated as described in the last section. 
In addition to the SPH calculations we have plotted in Figure [7] an estimate of the variation 
of the stride length with viscosity based on an analytical result obtained by Taylor (1952) 
in his discussion of the swimming of long slender bodies. A curve was fitted for the stride 
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length of the form 



S 



gi 

,.1/2 



+ S 2 , 



(6.4) 



where si and s 2 are arbitrary constants (since we have fixed 0) determined by fitting to 



our data set. The form of (6.4) is determined from (7.12), which will be discussed in the 
next section. For the present case we determine the coefficients using two points from the 
SPH results. We find S\ = 0.05774 and s 2 = —0.4386. The curve shows good agreement 



in the higher viscosity range where v > 5 x 10 4 and 3? < 500. We do not expect (6.4) to 
be valid for very high Reynolds number. 
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Figure 7: Stride lengths, in terms of length parameter a, for the motion in a fluid with different viscosities 
v but constant gait and mass. The crosses denote the SPH results and the filled circle denotes the result 
of Kanso et al. The dashed curve is based on an analytical result by Taylor, which is not expected to be 
valid for very high Reynolds number. Note that if v = 1CP 3 then Re = 250. 
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7 Speed and power scaling relations 



The speed with which the linked bodies move through the fluid depends upon a number 
of parameters. As already seen, the speed depends (at least) upon the mass of the bodies 
and the fluid viscosity. Additionally, we expect the speed to depend upon the ratios a/b, 
a/c, the frequency u>, and the amplitude (3. Similarly, we expect the power expended to be 
dependent on these parameters. Because the fluid is treated as slightly compressible there 
is a further non dimensional quantity ua/c s typically equal to 1/20 in our calculations. 
We neglect contributions from this quantity. 

The speed of the linked bodies, V, was estimated from the stride length divided by 
the time taken to complete the stride. Following the approach from the previous section, 
we take the stride length to be the average of the second and third strides. The time 
taken to complete a stride is 2tt/uj. 

In a biological creature the power expended for locomotion is provided by the ac- 
tions of the muscles which themselves depend on their body chemistry. In the case of a 
marine robotic vehicle the energy is provided by the engines within the vehicle. In our 
formulations the energy can be estimated from the constraint forces in the equations of 
motion. We calculate the average power V expended by the linked bodies over the time 
interval t± to t 2 from the expression 



V =J— I \ > ,V,-F,+ > />,-,. 1 dt. (7.r;) 
where and are the constraint forces and torques on body k respectively. Substituting 



the constraint forces and torques (2.10 2.15) into (7.5) gives 

v = — !— r (\ { e\n 2 - n x ) + xf\n 3 - n 2 j) dt. (7.6) 

£2 — h J tl v / 
V was calculated by numerical integration over the time of interest (in this case, the 

time taken to complete the first three strides). 

To simplify the velocity and power relations, we study the motion with bodies of 

fixed mass, fixed lengths a, b, and c and a fixed periodic-domain size. We then expect the 

speed to be given by an expression of the following form 

V = uaF^{3). (7.7) 
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Since the work done by the constraints is proportional to u 2 a 2 , and t oc 1/cu, we expect 
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Figure 8: Velocity of the linked bodies for different frequencies with constant Reynolds number and 
amplitude. The SPH results are shown by the crosses and the continuous curve is a best fit straight line. 

a power relation of the following form 

V = (3 2 u 3 a 2 G{^,(3). (7.8) 

It is useful to compare these scaling relations with the analysis of Taylor (1952). He 
shows that an infinite, flexible cylinder, along which a wave of amplitude B and wavelength 
A propagates with speed U = \u/(2n), moves with an average forward velocity V in a 
viscous fluid given by 

CnR[ /2 C(a) = 5 ffi a ffi } - J^W 2 ^ Ah{a) + 4Js(a)) ' (7 ' 9) 

where n = V/U, Co is a drag coefficient, and the functions 7, C, ii, I2, and ^3 are given 
by integrals. The quantity a is given by tana = 2nB/\ and when a is small enough 
a ~ 2-kB/X. For our oscillating bodies A ~ 6a, so that a ~ B/a. Because B/a is close 
to our amplitude (5 we replace a by (3 to convert Taylor's formula to a form appropriate 
for our system. Ri = Ud/v is a Reynolds number where the characteristic length is d 
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the diameter of the cylinder. Taylor's formula is an example of the relation (7.7) with a 
replaced by A. Our oscillating ellipses are similar to a small section of Taylor's oscillating 




Figure 9: Power expended by the linked bodies for different frequencies with constant Reynolds number 
and amplitude. The SPH results are shown by crosses and the continuous curve is a best fit cubic. 



cylinder and this suggests that Taylor's formula might provide a useful model for the 
scaling relations appropriate to the linked ellipses even though his calculations are for 
motion in three dimensions and ours are for motion in two dimensions. To that end we 
replace R\ by our Reynolds number 3? and expand Taylor's formula assuming n is small. 
We then find, with U oc ua, that 

— = A(aW 1/2 + B(a) (7.10) 

where the left hand side is the stride length scaled in units of a and a constant of pro- 
portionality has been absorbed into A and B. In the following we replace a by j3. The 
continuous curve in Figure [7j which applies to the case of constant (3 and u, was obtained 
by fitting the parameters A and B using two values of the viscosity. It can be seen that 
this gives a good fit to our results except at the lowest values of the viscosity coefficient. 
Furthermore, for constant viscosity and amplitude, the variation of V with frequency 
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deduced from (7.10) is 



V = c l u(uj l l 2 + c 2 ) 1 (7.11) 

where C\ and c 2 are constant for fixed amplitude, viscosity and length a. 

If a <C 1, the integrals in Taylor's formula can be expanded in a power series in a. 
If this is done, and we replace a by (3, we find that 

V 



£^1/2^5/2 + ^2 



(7.12) 



where k\ and k 2 are constants. These formula, with suitable values for the constants, give 
a very satisfactory fit to our results. 

Taylor (1952) also derives an expression for the power generated. If n <C 1, and the 
amplitude (3 < 1 we can expand Taylor's formula to find (we replace his notation W for 
the power by V) that 

^ = ^$(^ 1/2 P 2 + ^ 5/2 )- (7-13) 
When the viscosity and the lengths a, b and c are constant we can write this as 

v = c , yp 2 + c , y/ 2 p 5 / 2 . 



(7.14) 



Taylor's expressions for the velocity and power agree, as expected, with the general scaling 



relations (7.7) and (7.8). 



7.1 Speed and power for constant Reynolds number 

From the scaling relations we expect for a fixed amplitude and Reynolds number that 
V oc u and V oc u 3 . In order to test these relations the mass of each ellipse was kept 
constant (neutrally buoyant), with constant lengths a = 0.25, b = 0.2a and c = 0.2a and 
Reynolds number 200, with (3=1. The frequencies were in the range 0.5 < u < 2.0. 
In order to keep the Reynolds number fixed, the fluid viscosity was changed with the 
frequency of oscillation according to v = 1.25 x 10~ 3 u;. Additionally, the speed of sound 
required for the equation of state was constant with the value c s = 10 x 2a x co> max = 10 
where co> max = 2 is the maximum frequency used. 

Figure [8] shows that there is an approximate linear relationship between the velocity 



and the frequency, which is in substantial agreement with (7.7) except that the velocity 
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vanishes for u ~ 0.2. Figure M shows the power V against uj. The continuous curve is a 



cubic polynomial in agreement with (7.8). 



7.2 Speed and power for constant viscosity 

For the case where the viscosity is fixed we can be guided by Taylor's formula. We ran the 
calculations, again with neutrally buoyant bodies, and the same body length parameters. 
We chose v = 1.25 x 10~ 3 , such that 9ft = 200co>. We ran the calculations for a number 
of frequencies and amplitudes in the range 0.5 < u < 2.0 and 0.3 < (3 < 1.6. Note that 
for (3 > tt/2, the angle between two consecutive bodies is less than 90°. This case has no 
physical analogue to the long slender body of Taylor. 



The speed is plotted against u for three amplitudes 0.5, 0.9, and 1.3 in Figure 10 



We used two data points from the (3 = 0.9 data set in order to determine the constants 



k\ and k 2 in (7.12) which can then be written. 



V = 0.043396^ 3/2 /3 5/2 - 0.013039w/3 2 . 



(7.15) 



As Figure 10 shows, this gives a good fit to the (3 = 0.5 and 0.9 data. The (dashed) curve 



for (3 = 1.3 does not agree with the SPH results, but we do not necessarily expect it to, 



since (7.12) is only valid for (3 < 1. The velocity is plotted against (3 for fixed frequencies 



in Figure 11 The curves on these plots are again from (7.12), using the same values for 



k\ and k 2 as determined previously. Once again we see that the curves give a good fit for 
(3 < 1. We found that the peak velocity is achieved with (3 ~ 1.4 for all of these cases. 
The velocity appears to be smaller when the bodies swing through 90°, or more, relative 
to one another. 



The average power is plotted against u for fixed amplitudes in Figure [12} As with 
the velocity, we chose two data points from the j3 = 0.9 data set in order to determine 



the constants c" = 4.469 and c 2 = 0.798 in (7.14). This gives a very good fit to the SPH 



results. The power is plotted against (3 for fixed frequencies in Figure 13 Again, (7.14) 



with the same values for d[ and c' 2 ' gives a very good fit to the SPH results. 
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Figure 10: Velocity curves of constant viscosity and amplitude against frequency. Crosses, open circles 



and filled circles are for (5 = 0.5, 0.9 and 1.3 respectively. The curves are from (7.121 with constants 
fitted from the set with f3 — 0.9. The dashed curve is for (3 = 1.3, which is outside the range for which 



(7.121 is valid. 
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Figure 11: Velocity curves of constant viscosity and frequency against amplitude. Crosses, open circles 



and filled circles are for ui = 0.75, 1.25 and 2.0 respectively. The curves are from (7.121, which is 
appropriate for (3 < 1. 
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Figure 12: Average power curves of constant viscosity and amplitude against frequency. Crosses, open 



circles and filled circles are for f3 = 0.5, 0.9 and 1.3 respectively. The curves are from (7.141 




Figure 13: Average power curves of constant viscosity and frequency against amplitude. Crosses, open 



circles and filled circles are for ui — 0.75, 1.25 and 2.0 respectively. The curves are from (7.141 
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7.3 An application to the swimming of a leech 

Taylor applied his formula to estimate the speed of a leech which swims with shape changes 
shown in Figure 14 These shape changes only roughly approximate a sine wave down a 
very long cylinder since Taylor does not include end effects and, in any case the motion is 
only approximately sinusoidal (see for example frame 4). We can estimate the speed by 
making appropriate adjustments to our scaling formulae. Since these mimic Taylor's the 
only issue is whether the constants calculated from our simulations allow us to calculate 
the speed for an animal with different length scales. 

Returning to the relation (7.12), and determining the constants by fitting to our 
results for we find 

V = 0.01227 ^-j=) u 3/2 /3 5/2 - 0.05216c;a/5 2 . (7.16) 



The leech Taylor considered was approximately 0.08 m long, and travelled with a 
velocity 0.043 m/s. The gait of the leech produces a wave like motion along its body with 
an average speed 0.153 m/s which we estimate as being equivalent to 

U=^-\~uja, (7.17) 

where we take the wavelength as the length of the leech and set this to be approximately 
6a, the total length of our ellipses. This gives u ~ 12 s _1 . Taylor finds that the average 
value of B/X = 0.089 so that a = 0.56. Consistent with out earlier discussion, we 
can take (3 = 0.56. The thickness of the leech changes as it moves with an average 
diameter estimated by Taylor to be 0.0055m so that, for our ellipses, we can estimate 
b = 0.0027. The estimate of a for the leech is 0.08/6 and the ratio a/b ~ 5 as in 
our simulations. The viscosity coefficient is 10~ 6 m 2 /s. Substituting these values into 



(7.16) with (3 = 0.56, we find V = 0.037 m/s, which compares favourably with Taylor's 
estimate from the experiment of 0.043 m/s. It might be thought that this agreement is a 
lucky coincidence but, taken with the agreement of our results with formulae modelled on 
Taylor's expression, it does suggest that the speed of a three dimensional long thin body 
is similar in form to that of three linked, long, thin ellipses in two dimensions. 
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Figure 14: Frames showing the motion of a leech against a background of squares of side 2 cm. The time 
interval between frames is 15 s. The image is part of Figure 7 in Taylor's paper. The shape of the leech 
is roughly similar to that of our three oscillating ellipses. 



8 Optimal motion with constant viscosity 

We now ask the question: for a given fluid viscosity, body size and body mass configura- 
tion, what is the frequency and amplitude to move with a given speed, while expending 
the least power? We do not attempt to obtain the optimum frequency and amplitude for 
all possible ellipses. Instead we have the more modest aim of determining if an optimum 
set of parameters exists for a typical set of ellipses. With this in mind we consider neu- 
trally buoyant bodies, with a = 0.25, b = 0.2a and c = 0.2a. The kinematic viscosity 
is v — 1.25 x 10~ 3 , so that 9ft = 200cj. In Figure 15 we show the contours of constant 



velocity and constant power in the (u, f3) plane. These contours were obtained by using 
the results of the simulations to fit the velocity with polynomials of the form 

2 4 

i=i j=i 

and the power V with the function cf3 2 uj z . While (8.1) does not include the fractional 
powers we derived by comparison with Taylor's formula it is still possible to a satisfactory 
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over the whole range of uj and (3. It is clear from Figure 15 that there is a set of values of 
uj and (3 which will give a specified speed with minimum power. The dashed line in Figure 



15 gives the optimum set of uj and f3. This line was calculated by traversing a contour of 



constant power and finding the uj and f3 which give the maximum speed. It is interesting 
to note that the optimal motion is close to (3 ~ 1.2 regardless of the frequency. 

We can estimate some properties of the contours in Figure [15] without detailed 



numerical calculations. For constant viscosity, we can estimate from (7.12), for fixed 
lengths a, b and c that 



/3 2 oc 



V 



(8.2) 



uj^uj 1 / 2 + k[) ' 

where we have replaced /3 5//2 by (3 2 which is reasonably accurate for 0.5 < f3 < 1, and k[ 



is a constant. Similarly from (7.14) we estimate 

(3 2 oc V 



(8.3) 



W 5/2( w l/2 + jfc /)- 

Where k' 2 is a constant. These expressions show that (3 increases faster as uj decreases for 
the constant V curves than for the constant V curves when 0.5 < (3 < 1. This gives the 



shape of the contours on the left hand side of Figure 15 



9 Conclusions 

The principal results of this paper are (a) that the accuracy of the SPH algorithm for 
linked bodies moving in a fluid has been established, (b) that the variation of the calculated 
speed and power output take simple forms consistent with scaling relations, (c) that there 
is remarkable agreement between the two dimensional results and those Taylor obtained 
for the swimming of long narrow animals in three dimensions, and (d) the minimum power 
to produce a specified speed for a given gait has been calculated and forms a basis for 
other such calculations. 

The first of these results has been obtained by resolution studies and by comparison 
with the results of Eldredge (obtained for massless bodies in a viscous fluid) and Kanso 
et al. (obtained for neutrally buoyant bodies in inviscid fluids). In both cases the relevant 
results are limits of our calculations. In the case of Eldredge we estimated his value from 
a series of calculations where the mass of the body was changed. In the case of Kanso et 
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Figure 15: Velocity and average power contours on a frequency- amplitude plot. The solid lines are the 
velocity contours, and the dashed lines are the power contours. The thick dotted line is the curve of 
optimal motion. 
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al. the viscosity was steadily decreased so that the Reynolds number increased from 50 
to 5000. 

The second and third results were obtained by fixing the dimensions of the bodies 
and their masses but changing the frequency and amplitude of the gait. The simulations 
show that the the results have a simple dependence on frequency and amplitude which is 
similar to that found by Taylor (1952). These results suggest that the drag forces on a 
long thin ellipse in two dimensions is similar to that on a cylinder in three dimensions. 
In particular, it suggests that the drag has two additive contributions. One varying with 
Reynolds number as l/K 1 / 2 and one depending on the square of the velocity relative to 
the fluid. We are unaware of calculations or analysis which would confirm this conjecture 
in detail. It is clear however, that there will be pressure forces proportional to the square 
of the velocity on the bodies, and viscous forces due to flow along and between the ellipses. 

The result (d) shows that the efficiency is poor if the linked bodies are driven with 
a gait amplitude which is too large or too small. We find that the optimum performance 
occurs when the amplitude (5 ~ 1.2 or, equivalently, when the angles between the links 
varies between ±7r/3. 

The formulation we have used is general and can be immediately applied to the 
motion of linked bodies in stratified fluids, or with a free surface, or within complex 
boundaries, or with more complex constraints including those where the gait depends on 
the positions of the bodies in the domain. We are currently studying these problems. 
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